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Summary 

Incompressible two dimensional calculations are reported for the impulsively started 
lid driven cavity with aspect ratio two. The algorithm is based on the time dependent 
streamfunction equation, with a Crank-Nicolson differencing scheme for the diffusion terms, 
and with an Adams-Bashforth scheme for the convection terms. A multigrid method is 
used to solve the linear implicit equations at each time step. Periodic asymptotic solutions 
have been found for Re = 10000 and for Re = 5000. The Re = 5000 results are validated 
by grid refinement calculations. The solutions are shown to be precisely periodic, and care 
is taken to demonstrate that asymptotic states have been reached. A discussion is included 
about the indicators that are used to show that an asymptotic state has been reached, and 
to show that the asymptotic state is indeed periodic. 


*Work funded under Space Act Agreement C99066G. 



1: Introduction 


Cavity flows have been a subject of study for some time. These flows have been 
widely used as test cases for validating incompressible fluid dynamics algorithms. The 
preponderance of such studies have addressed the steady flow problem, such as Ghia et. 
al. [4] and Schrieber and Keller [17], and the potentially rich unsteady cavity dynamics 
have only recently begun to be addressed. Greatly increased computational capabilities 
make possible the study of unsteady flows and initiates a new chapter in numerical analysis, 
determining the qualitative properties of the solutions of time dependent partial differential 
equations from their simulation. Examples of this type include a study of transition to 
turbulence in a two dimensional cascade flow by Fortin et. al. [3], and a study of the 
transition to chaos for a Boussinesq fluid in a vertical cavity by Paolucci and Chenoweth 
[16]. Recent work on the numerical analysis of large time and asymptotic solutions of 
partial differential equations includes the study of the convergence of attractors for finite 
dimensional approximate systems to the attractor of the original system by Hale et. al. 
[11], the study of the large time behavior of Galerkin approximations to Navier-Stokes 
equations by Constantin et. al. [2], and the study of finite element approximations of the 
nonstationary Navier-Stokes equations by Heywood and Rannacher [12]-[14], and especially 
for behavior as t — > oo in [14]. A theoretical treatment of the dynamical systems approach 
to the Navier-Stokes equations can be seen in Constantin and Foias [l], and in Temam 
[19]. Both the qualitative nature of asymptotic states, and the way in which they develop 
are new and exciting areas of research in partial differential equations. The driven cavity 
problem gives a computationally reasonable model fluid problem in which to investigate 
the qualitative features of the solution space of a physically reasonable infinite dimensional 
dissipative dynamical system. In contrast with other unsteady problems that are currently 
being investigated, the driven cavity is a self contained system with realistic nonperiodic 
boundary conditions, with steady nonperiodic forcing, and without artifical throughflow 
boundary conditions. 

The vortex dynamics for unsteady Navier-Stokes flow in the driven cavity for various 
aspect ratios and Reynolds numbers with an impulsively started steady lid were studied in 
Gustafson and Halasi [8] and [9], and in Goodrich and Soh [7]. A persistant final oscillation 
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in the aspect ratio two driven cavity with a 40 x 80 grid was shown in [9] at Re = 10000 
and t « 300, and the conjecture was raised that a qualitative transition had occurred to 
an unsteady asymptotic flow. The present paper continues that study by addressing the 
question of whether an unsteady asymptotic solution exists in the driven cavity with aspect 
ratio two. In section 2 we discuss some indicators for monitoring the development of a flow, 
for identifying an asymptotic flow state, and for presenting the qualitative nature of that 
state. In section 3 we present and discuss the algorithm that is used for the computations 
reported in this paper. In section 4 we discuss a periodic solution at Re = 10000 on a coarse 
48 X 96 mesh, comparing with and extending the results in [9]. In section 5 we present 
and discuss the main results of this paper, two periodic solutions at Re = 5000 on 48 X 96 
and 96 x 192 grids, with periods of approximately 2.477 and 2.309, respectively. The two 
asymptotic states computed at Re — 5000 have qualitatively similar dynamics. All three 
of the computed asymptotic states reported in this paper are periodic with very precise 
repetition of the asymptotic cycles. Based upon these computations we conclude that a 
Hopf bifurcation does occur in 
number Re c < 5000. In section 


the aspect ratio two driven cavity for a critical Reynolds 
6 we discuss a number of continuing and remaining issues. 


3 



2: Measures of Qualitative Flow Features 

Flow visualization is a fundamental problem for computational fluid dynamic is ts as 
well as for experimentalists. This is especially true for the resolution of small scale flow 
details. A primary issue for the simulation of flow evolution to an asymptotic state is to 
identify measures that will give reasonable assurance that an asymptotic state has in fact 
been attained. This issue is somewhat obscured when investigating flows that converge 
to a steady state. We will list and briefly discuss the measures that we have used for the 
qualtitative representation of flow simulations. 

2.1: Field Representations 

The first series of measures and indicators represent the standard two dimensional 
flow field data. The representations that we have found useful are: 

1] Streamfunction Contour Plots - These are useful for defining the large scale flow 
features. When overlayed with normalized velocity vectors, they give a good balance 
between large and small scale resolution. 

2] Streamfunction Surface Plots - These plots are good for presenting large scale flow 
features. Like all surface plots they give an immediate and vivid impression of scaling 
throughout the flow field, but depending upon the view of the surface, prominant flow 
features can mask smaller scale structures. 

3] Velocity Vector Plot - These plots are standard and are almost indispensible for ac- 
tually visualizing the flow. Normalized velocity vectors are better for revealing the 
smaller scale structures. 

4] Kinetic Energy Contour Plot - These show the momentum scale for the entire flow. 

5] Kinetic Energy Surface Plot - These are extremely useful for showing the momentum 
scale for the flow and its various features. This is an excellant complement to either 
a normalized velocity vector plot or a streamfunction contour plot overlayed with 
normalized velocity vectors. 

6] Vorticity Contour Plot - These are extremely useful for following vortex dynamics, 
particularly for vortex pairing and splitting as in [10]. Vorticity can be concentrated, 
so vorticity plots often reveal less than streamfunction contours about overall flow 
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dynamics. Vorticity data is very useful for highlighting the resolution failures of an 
algorithm. 

7] Vorticity Surface Plot - Since vorticity contours tend to be concentrated, a surface 
plot can give a better global view of the vorticity for a flow. Large regions of constant 
vorticity contrast well with local singularities and steep gradients. 

8] Pressure Gradient Vector Plot - These clearly reveal principle vortex centers and their 
relative significance as field sources compared to secondary vortices. 

The algorithm that we use for the results reported in this paper does not require a pressure 
solution, so we will not present any pressure data. The pressure gradient has been used to 
good effect by Gustafson and Halasi in [8]. All of the field representations that do not use 
pressure data were produced for this study, but for the sake of brevity the flow field will 
be represented just by streamfunction contour plots with and without normalized velocity 
vector overlays, and by surface plots of streamfunction, kinetic energy, and vorticity. 

2.2: Indicators of Dynamics 

The second series of measures and indicators are for tracking convergence to an asymp- 
totic state, and for understanding the qualitative nature of that state. These various 
measures and indicators are scalars that are either global indicators or point data values. 
Global indicators can be obtained by using various mathematical norms on the flow data, 
or on the flow change data. Any of the scalar valued data can be obtained at each time 
step, and there are a variety of ways to display the information content of such scalar time 
series data. A simple graph of data plotted against time gives a quick and familiar view 
of the time changes in scalar data, although this type of presentation can mask important 
details. A spectral density or power spectrum is a well known device from signal processing 
which identifies the flow frequencies and their relative strengths. A phase diagram can be 
presented for scalar data either for one variable with the value at time t n plotted against 
the value lagged by k time iterations at , or for two separate variables at the same time 
step plotted against each other. Spectral densities and phase diagrams are widely used 
to investigate and present the qualitative features of nonlinear dynamical systems. Both 
of these methods are particularly useful for classifying unsteady flows as either periodic 
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or aperiodic, for identifying critical parameter values for transition between qualitatively 
different flows, and for identifying chaotic or turbulent flows. The indicators that we have 
found useful are: 


9] The Relative Li Norm of the Streamfunction Change Per Time Step - This is a global 
measure, and is calculated as 


E, ,, l <; 1 - Kt 

e„i*ti 


This measure is similar to standard convergence indicators, and is easy to compute, 
but it is an integrated measure and tends to mask relatively small scale flow features. 
We did not look at the combined relative change in streamfunction and vorticity, which 
could be an interesting and more sensitive measure than just the streamfunction or 
vorticity data considered separately. 

10] The Maximum and Minimum Streamfunction Value - These measures provide useful 
information about the location of the two main vortices in the cavity, and about the 
dynamics of the asymptotic solution at the center of each of these vortices. In addition 
to the global extremes of the streamfunction, it would be interesting to also record 
other local extremes in order to identify and track vortex centers, and to give their 
relative intensities. 

11] The Relative L x Norm of the Vector Field Change Per Time Step - This is also a 
global convergence measure, and is calculated as 


E,. i (Kr-«r.,l + Kr-<,l) 

E,.,(Kn+Kri) 


We used this measure along with the relative Li norm of the streamfunction change 
data, and found that they behaved very similarly to each other. A choice between 
these two measures can be easily made on the basis of which variable set is being used. 

12] The Total Kinetic Energy - Calculated for the entire grid minus the boundaries at 
each time step as 

^AxAvY, K, + 1 C- 
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This measure was more sensitive than the streamfunction change data and was a 
very useful complement to that data. The time dependent dynamics of a flow are 
formulated in terms of momentum balance equations, so that such a global momentum 
measure would seem to be very appropriate data for the assessment of the evolution 
of a flow. In fact, our experience (see section 5) has indicated that very significant 
evolution can still occur even though the streamfunction or velocity time step change 
data becomes as small as 0(lO~ 7 ) or 0(1O” S ), with the total kinetic energy remaining 
noticably sensitive to the continuing evolution of the flow. This experience suggests 
that co mm on practices for stopping calculations could produce misleading results 
about the qualitative features of asymptotic flows. 

13] The Maximum Acceleration - Determined for the entire interior grid as 


u 


1 


u; 


max< 




At 


M- 


The point where the maximum acceleration occurs can change with the time step, 
possibly in a very discountinuous way. Consequently, this indicator can be sensitive 
to the time scales and dynamics of both the local change per time step in the point 
values of the velocity field, and the global convective transport processes within the 
flow. There can be subtle and complex relationships between these processes. Note 
that this indicator is not the maximum convective derivative, which would be easy to 
obtain and could be interesting to try. 

14] Streamfunction Value at a Point 

15] Velocity Component at a Point 

16] Kinetic Energy Value at a Point 

17] Vorticity Value at a Point 

For the calculations reported in this paper we found the relative L x streamfunction change 
norm and the total kinetic energy to be very good complements for indicating convergence 
to an asymptotic state. The streamfunction change per time step data showed convergence 
earlier, and the total kinetic energy showed final convergence to the asymptotic state. This 
conclusion was augmented by point data time series for streamfunction, kinetic energy, and 
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vorticity. The sensitivity of the indicators increased from the streamfunction, to the veloc- 
ity components defined as streamfunction derivatives, then to the kinetic energy defined 
with squared velocity components, and finally to the vorticity defined as the Laplacian of 
the streamfunction. 
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3: An Algorithm For Time Dependent Incompressible Navier-Stokes Equations 

This section briefly presents the implicit finite difference streamfunction algorithm 
for the unsteady incompressible Navier-Stokes computations reported this paper. Further 
specific algorithm details may be found in Goodrich [6] or Goodrich and Soh [7]. An explicit 
MAC (marker and cell) primitive variable scheme was used in [8] and [9], The time step 
restriction for such an explicit scheme is prohibitive for truly long time studies. The 
qualitative flow dynamics of both codes agree for the computations that are comparable. 

3.1: Flow Equations 

In a bounded open region fl, the general dimensionless Navier-Stokes equations for 
incompressible flows are 

— + (u • V)u — A u = - Vp + F, for x in 0, and t > 0, 

dt v ; Re 

V • ii = 0, for x in fl, and t > 0, 


u(x,0) = a(x), for x in fl, at t = 0, 

i?[u(x,f)] = b(x,t), for x in 3fl, and t > 0, 

where u is the velocity, p is the scaled pressure, Re is the Reynolds number, F is the 
volume force per unit mass, c?fl is the boundary of fl, and B is the operator that defines 
the boundary conditions. We will ignore the body force F. For the simulations of two 
dimensional cavity flows that we will be reporting, the computational domain fl is just a 
rectangle as in Figure 1. The length and velocity scales for these cavity flows are taken to 
be the lid length and velocity, so that the lid velocity is always 1. The velocity field for 
two dimensional flows may be written in terms of the streamfunction ifr as 

u(x,t) = — and v(x,t) = for x in fi and t > 0, (la) 

ay ox 

and the dimensionless equations for evolution of time dependent viscous incompressible 
flows may be written as the streamfunction equation 


d A 0 
dt 



— A — , for x in fl, and t > 0. 
dydx 


( 16 ) 
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In this formulation the data for the impulsively started driven cavity consists of the initial 
values 

V>(x,0) = 0, for x in 0, at t = 0, (lc) 

and the standard boundary conditions (see Figure 1) 


V>(x, t) = 0, for x in 3H, and t > 0, 


(Id) 


dtp 

drj 

dtp 


(x,t) = 1, for x on the cavity lid, and t > 0, 


M 


— -(x,t) = 0, for x on the cavity walls, and t > 0, 

drj 

where •£- is differentiation in the exterior normal direction at the boundary. 
#*» 


3.2: The Discretization 


The discrete approximation for the streamfunction ip = {ip{x,t) : x € fi,f > 0} 
will be taken to be z = {i m : m = 0,1,...} = {«,”• : * = 0,1 J = 0,1,... J, m = 
0, 1, . . . } on the rectangular computational mesh. We will use a uniform grid on the 
rectangular domain, and we will use centered spacial differencing for equation (lb), with 
a Crank-Nicoison time differencing for the diffusion terms, and a second order Adams- 
Bashforth time differencing for the convection terms. Let La be the conventional five point 
centered difference approximation to the Laplacian, let Bi be the conventional thirteen 
point centered difference approximation to the Biharmonic operator, and let 6 X and 6 V be 
the conventional centered difference operators 






m 

Z i- hi 


, and )t\j 


y-n 


2Ax ’ wv “ ,,,J 2A y 

With this notation and differencing scheme, we may discretize equation (lb) as 


La(i " +,) ~ ^ Bi(£ " t,) 

= La(£”) + ^Bi(i") - ^14 


I. (<,(£-)La(£")) (<.(£" )La(£”)) 


( 2 ) 


+ Y )La(*-'))-<i (<.(£" -*)La(i-‘)) . 
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In order to incorporate the boundary condition (le), the standard finite difference ap- 
proximations are used for modifying the discrete equation (2) at grid points next to the 
boundary. Note that this algorithm is second order accurate in both time and space. The 
velocity components (u t "V , v ™ } ) are directly recovered from the discrete streamfunction 
solution as 


;«Vn - and. ~ z ?~ i.y)« 


Notice that the velocity components are both defined at each grid point and not at different 
locations on a staggered grid. In our formulation we use a grid centered central difference 
expression for the mass conservation equation. The velocity solution is exactly discretely 
divergence free with respect to this divergence operator. 


3.3: A Solver 


The operators in the implicit time stepping equation (2) do not change with time. In 
fact, the discrete streamfunction solution i n + l from equation (2) depends upon the time 
step only through the solutions at the previous two time steps. The problem that must be 
solved at each time step may be written as 


La(i ” +1) " ^Bi(£" +1 ) = f(i- f i-‘), 


where f is the discrete source term from the right hand side of (2). This problem is a 
discretization of 

At , 

A for x in fi, att = f„ + 1 , 


an elliptic partial differential equation that combines the Laplace and Biharmonic opera- 
tors. At each time step we use a multigrid method to solve these equations, as in Goodrich 
[6], The Biharmonic operator is factored as two Laplacians as in Linden [15], point Gauss- 
Seidel relaxation is used for the smoothing operator, and linear restriction and prolongation 
operators are used. A V-cycle iteration scheme is used, with 3 iterations per grid level 
while coarsening, and none while refining. At each time step, 10 to 15 iteration cycles are 
used to reduce the residuals in (3) to less than 5.0 x 10~ 12 . 
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4: Numerical Results At Re = 10000 

The first n um erical results that suggest a periodic solution for the driven cavity with 
aspect ratio two are in Gustafson and Halasi [9]. These computations are at Re — 10000 
for 0 < t < 360, using an explicit primitive variable Euler-MAC scheme on a 40 x 80 grid 
with At = The published results are a series of normalized velocity vector plots of 

the flow field that suggest an eventual periodic pattern of evolution and interaction in the 
cavity vortex dyn am ics. The grid resolution that was used at this Reynolds number leaves 
open the possibility that the computational mesh might be filtering out dynamical effects 
with smaller spatial scales. The reliance upon just the normalized vector plots leaves open 
the question of whether a periodic, a quasi-periodic, or an aperiodic flow was actually 
observed. The relatively short time interval that was calculated and the small number 
of flow periods that were observed does not ensure that the initial transients were fully 
dissipated and that a true asymptotic solution was obtained. 

An independent series of calculations was initiated to verify the results reported by 
Gustafson and Halasi [9]. These computations are at Re = 10000 for 0 < t < 1800, using 
the algorithm described in section 2 on a 48 x 96 grid with A t = . These computations 

show that the asymptotic state on this grid is a periodic flow, although the asymptotic state 
occurs only at a much later time than suggested by Gustafson and Halasi [9]. The results 
of these computations are similar to the data for the stronger results that we have obtained 
at Re = 5000, which are the focus of this paper, so for the sake of brevity the data at 
Re = 10000 will be discussed but will not be shown. Several indicators were used to ensure 
that an asymptotic state was actually reached. The first global indicator is the relative L x 
norm of the streamfunction change per time step. This indicator suggests that a periodic 
asymptotic state has been reached at t « 300, as reported in Gustafson and Halasi [9]. 
The global maximum and minimum streamfunction values were tracked as the second and 
third global indicators of approach to an asymptotic state. The observed data showed that 
as the flow evolved, the maximum and minimum streamfunction values on the grid begin 
to oscillate at the centers of the upper and lower vortices. The maximum streamfunction 
value in the secondary lower vortex reaches its asymptotic oscillation at t « 550, and 
the minimum streamfunction value in the primary upper vortex reaches its asymptotic 
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oscillation at t « 900. A fourth global indicator is the relative L t norm of the vector field 
change per time step. This indicator suggested that a periodic asymptotic state had been 
reached at t « 500. The fifth global indicator is the total kinetic energy for the entire grid 
minus the lid. This data reaches a periodic asymptotic state at t « 1500. Besides this 
global data, local flow data was gathered at selected fixed points in the grid. At points 
with coordinates (t, j) = (44, 50) and (t,j) = (16,46) the values of the streamfunction, the 
x and y velocity components, and the vorticity were recorded, and the kinetic energy at 
the point was calculated, at each time step for 1600 < t < 1800. The first point is near the 
end of the wall jet that descends from the lid along the upper half of the downstream wall. 
The second point is near the end of the bounded shear layer as it approaches the wall after 
crossing the middle strip of the cavity. This data shows clean spectral signatures with one 
fundamental frequency. The data at the end of the bounded shear layer shows energization 
of harmonics of the one fundamental frequency as a result of convection through the shear 
layer. Based upon all of these indicators the discrete dynamics on this grid do appear 
to reach a periodic asymptotic state at t « 1500, with a period between 223 and 224 
time steps, or between 3.48 and 3.50 nondimensional time units. The periodic asymptotic 
flow on this 48 x 96 grid at Re = 10000 is qualitatively the same as the vortex dynamics 
observed in [9]. 

A natural question is whether or not the periodic flow for this discrete system is a 
grid dependent phenomenon. This is a particularly relevant question for this simulation 
since the best steady solutions for the square cavity at Re = 10000 are on a 256 x 256 
uniform grid, and the 48 x 96 uniform grid that we are using clearly does not give anywhere 
near the same resolution. Another serious cause for concern is that the results reported by 
Gustafson and Halasi [9] on a 40 x 80 grid suggested a period of approximately 4.5, while 
our 48 X 96 grid calculation gives a period between approximately 3.48 and 3.50. To begin 
checking for grid dependence we conducted a calculation at Re = 10000 for 0 < t < 500 
using the algorithm in section 3 on a 96 x 192 grid with A t = . This calculation 

shows a much more complex time evolution than the results for the same time interval 
on the 48 x 96 grid, and it is possible that the 96 x 192 grid simulation is itself masking 
even finer scale dynamics. This suggests that the coarser grid was filtering smaller scale 
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dynamic processes that could contribute to a more complex asymptotic state than a single 
frequency periodic flow. We have in fact concluded that the Re — 10000 simulation on 
a 48 x 96 grid is most likely inadequate for resolving all but the coarsest features of the 
cavity flow, and the periodic discrete solution on the coarse grid at Re = 10000 cannot be 
accepted as giving a reliable portrait of the actual asymptotic continuum flow dynamics. 
Consequently, we discontinued the calculations at Re = 10000 and decided to initiate 
discrete simulations at Re = 5000 where less grid resolution would be required to capture 
the continuum dynamics. 
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5: Numerical Results At Re = 5000 


Two different simulations were computed to asymptotic periodic states at Re = 5000 
for the impulsively started driven cavity with aspect ratio two. Both simulations use the 
algorithm described in section 3. The first simulation uses a 48 x 96 grid with A t = — 

64 5 

and the second simulation uses a 96 x 192 grid with At = The simulations at 

e 10000 and Re — 5000 are similar in the sense that they all demonstrate a very clean 
and qualitatively similar periodic solution to the nonlinear discrete dynamical system that 
is being used as an approximation to the fluid flow in the cavity. The conclusion of the 
discussion in Section 4 about the Re = 10000 simulations is that a 48 x 96 grid is inadequate 
for resolving all but the coarsest features of the cavity flow, so the calculated periodic 
discrete solution on the coarse grid at Re = 10000 cannot be accepted as giving a reliable 
portrait of the actual asymptotic continuum flow dynamics. The Re = 5000 simulations 
are essentially different in this regard, and we may reasonably conclude that they are 
likely to represent at least the qualitative dynamic features of the actual continuum flow 
dynamics. The first supporting evidence is that the refined 96 x 192 grid has a resolution 
close to the 128 x 128 grid used in Goodrich [6] with the same algorithm to produce a 
steady square cavity solution in good agreement with the standard published Re = 5000 
solutions. The second supporting evidence is that the grid refinement calculation produces 
flow dynamics that are close to the coarse grid results. In particular, the period given by 
the coarse grid simulation is 2.469 < T < 2.484, while the refined grid simulation gives 
2.305 <T< 2.313. 

5.1: Re = 5000 on a 48 x 96 Grid 

The computation at Re = 5000 on the coarser 48 x 96 grid will be presented with only 
a small amount of detail. The computation is for 0 < t < 4100, and the flow reaches a 
periodic asymptotic state at t ft 1100 with a period between ft 2.469 and ft 2.484. 
This flow is completely periodic with no apparent secular trend for 1100 < t < 4100, or 
for approximately 1200 oscillation cycles. 

The entire flow field at t = 4000 is shown in Figure 2 as a streamfunction contour 
plot. The main features of the flow are the primary and secondary circulations in the 
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upper and lower cavity, the wall jet that descends from the lid along the upper half of the 
downstream waU, and the wavey shear layer between these two main circulations. There 
are tertiary corner vortices in the two lower corners that are separated from each other by 
one sixth of the bottom wall where the secondary vortex in the lower cavity is attached 
to the lower boundary. Surface plots of streamfunction, kinetic energy, and vorticity show 
some lack of grid resolution especially in the boundaries near the end of the wall jet, and 
in the upper left corner. The flow at this instant is generally representative of the flow at 
any time during the entire cycle of the periodic asymptotic state. The vortex dynamics 
for this periodic asymptotic solution are similar to the dynamics for the refined 96 x 192 
grid. A key feature of the periodic vortex dynamics to the appearance of a pair of sma 
counterrotating tertiary vortices along the downstream wall slightly below the end of the 
wail jet descending from the lid. The stronger and higher tertiary vortex periodical y 
appears below the wall jet and is converted as a wavy disturbance along the shear layer 
between the primary and secondary vortices in the cavity. More detail will be presented 

for the similar dynamics of the refined grid solution. 

Figure 3a shows the spectral signature of the streamfunction at the point with coor- 
dinates (i,/) = (44,50) or (*,») « (0.92,1.04), and Figure 3b shows this data as a phase 
portrait plot of bo«* <000 < t < 4100. Similar spectral signatures and 

phase portraits are shown for the kinetic energy in Figures 3c-d, and for the vorticity in 
Figures 3e-f. This point is near the end of the wall jet descending from the lid along the 
downstream waU, and the data represents a sample with 6400 time points, with between 
158 and 159 data points in each period, and for more than 40 periods. Figures 
similar data recorded at a point near the end of the bounded shear layer across the mid- 
dle of the cavity with coordinates (.', j) = (16,46) or (*,„) « (0.33,0.96). The values of 
streamfunction, kinetic energy, and vorticity from these two points aU have extremely clean 
spectral signatures with one fundamental frequency. Harmonics are energised especially 
for vorticity by the effect of the bounded shear layer. The narrowness of the plot lines 
indicates the precise repetitiveness of this data for these 40* periods. This point data 
indicates a precisely periodic solution at Be = 5000 on the 48 X 96 grid. 
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5.2: Re = 5000 on a 96 x 192 Grid 


The computation at Re = 5000 with the finer 96 x 192 grid is the main result and 
will be presented in some detail. The computation was for 0 < t < 4100, and the flow 
reaches a periodic asymptotic state at t fa 3700 with a period between « 2.305 and 

« 2.313. This flow is completely periodic with no secular trend for 3700 < t < 4100, 
or for more than 170 oscillation cycles. 

A plot of the relative L x norm of the streamfunction change per time step is shown 
for 5 < t < 100 in Figure 5a, for 100 < t < 4000 by intervals of 300 nondimensional time 
units in Figures 5b-n, and for 4000 < t < 4100 in Figure 5o. The total kinetic energy is 
shown for 5 < t < 100 in Figure 6a, shown by intervals of 300 nondimensional time units 
for 100 < t < 4000 in Figures 6b-n, and for 4000 < t < 4100 in Figure 6o. Note that 
the vertical scales in Figures 5 and 6 are all different, and that there are two time scales 
in each figure. The data in Figures 5 and 6 shows a dramatic history. Up to t ss 700 
the time history seems to be following a dissipating evolution. At t ta 700 the relative 
L x change per time step and the total kinetic energy both start to increase in absolute 
level with oscillations that increase in amplitude. There is also a beating process with the 
interplay of multiple frequencies. There is a peak in the global kinetic energy at t « 1600, 
and then the global kinetic energy decreases until t « 2600, with continuing increases in 
the amplitude of the kinetic energy oscillations. At t « 2600 the streamfunction time 
step change starts to decrease while the kinetic energy begins to increase again. Between 
t = 2700 and t — 2800 both variables start to oscillate in what seems to be an erratic 
manner until t « 2950. After this the data steadies down and becomes very regular until 
the streamfunction time step change data converges at t « 3250, and the kinetic energy 
data converges at t « 3700. Contour and surface plots for streamfunction, kinetic energy, 
and vorticity on the refined grid show no visually apparent changes for 500 < t < 2000. 
Noticable pulsating changes begin to occur in this surface data only for 2000 < t < 2500, 
along with a noticable increase in the momentum in the lower half of the cavity. Note that 
the relative L x norm of the streamfunction change per time step is 0(5 x 10 — 7 ) for t m 700. 
This is below what might easily be taken as a small test value used to stop the calculation 
with the result declared to be a steady asymptotic state. The actual asymptotic state has 
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extremely regular oscillations in the relative Li norm of the streamfunction change per 
time step, with relative L x norm values that are 0(3 x 10" 4 ) for t > 3700. 

The precisely periodic nature of the asymptotic flow is shown in Figure 7a by the 
spectral density of the relative L x norm of the streamfunction change per time step, and 
in Figure 7b by the spectral density of the total kinetic energy, both for 4000 < t < 4100. 
Both spectral signatures show a single fundamental frequency near 0.433 « . Note that 

the streamfunction change data has one harmonic. Both of these spectral signatures are for 
global data, and reflect the global dynamics on the entire grid. The global maximum and 
minimum streamfunction values were also recorded, and are shown together in Figure 7c 
as a phase portrait plotting ^ or 4000 < t < 4100. This plot is for data taken 

at 12800 discrete time steps, with between 295 and 296 time steps per cycle, for more than 
43 complete cycles. All of this data is plotted sequentially with time, so that the physical 
plot lines have been drawn over 43 times. The narrowness of the lines substantiates the 
precisely periodic nature of the asymptotic solution. The global minimum and maximum 
always occured at the centers of the primary and secondary vortices in the upper and lower 
halves of the cavity. The effect of the periodic asymptotic flow is still apparent even in 
the centers of the two largest vortices. The entire streamfunction surface is vibrating with 
extreme precision. Note that the two extremes are oscillating approximately one half of a 
period out of phase with each other. 

A qualitative portrait of the overall solution field at the single instant t = 4100 
is given in Figures 8a-d. The streamfunction surface plot in Figure 8a shows a wavey 
disturbance or swelling that is periodically convected by the bounded shear layer across 
the middle of the cavity. This particular view shows the disturbance just before it is 
absorbed and dissipated by the boundary layer flow at the end of the bounded shear layer. 
The kinetic energy surface plot in Figure 8b shows a periodic wave in the momentum 
ridge near the center foreground. The dramatic low point in the kinetic energy ridge 
is in the bounded shear layer just above the front of the streamfunction swelling, with 
local max ima and minima in the kinetic energy just upstream from this low point. There 
seems to be a periodic pulse of momentum associated with the periodic streamfunction 
disturbance. The vorticity surface is shown in Figure 8c. Note that the surface plots of 
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streamfunction, kinetic energy, and vorticity all show smooth resolution of the flow solution 
over the entire grid. The streamfunction contour plot in Figure 8d shows a recirculating 
boundary layer at the bottom of the cavity. This feature appears in the evolution of the 
flow before t = 500. Note that the streamline defining this feature is for ip = —0.001, 
while the streamline closest to the center of the primary circulation is for ip = —0.090. 
This recirculating boundary layer in the bottom of the cavity is imperceptible in the 
surface plots of streamfunction, kinetic energy, and vorticity, but it may be highlighted 
by a low enough streamfunction contour or a normalized velocity vector plot. There are 
two local weak tertiary recirculations diagonally out of the corners, within the ends of 
the figure eight shaped streamline for ip = —0.090, and above the recirculating boundary 
layer along the lower wall. It appears that a single large third vortex is partially formed 
in the bottom of the cavity, but that there is not quite enough room below the secondary 
vortex to allow for the joining of the two weak tertiary vortices above the lower wall. The 
refined grid calculation shows quartiary corner vortices, resolved with between one and four 
grid points. The recirculating boundary layer along the lower wall is not present in the 
coarse grid asymptotic flow, which has two tertiary vortices in the lower corners separated 
by a short stretch of grid points where the secondary vortex attaches to the lower wall. 
The asymptotic coarse grid solution has a smaller total amount of kinetic energy on the 
interior grid than the asymptotic flow on the refined grid. These differences between the 
two asymptotic flow solutions are probably caused by the better resolution of gradients 
and momentum diffusion on the refined grid. 

Figures 9a-k are a series of eleven plots at time intervals of 0.25 over slightly more than 
one complete periodic cycle starting at t = 4100. These eleven plots are streamfunction 
contours with normalized velocity vector overlays in a central band across the cavity, 
and they show the essential vortex dynamics in detail for a typical cycle of the periodic 
asymptotic flow. The first interesting feature of these dynamics is the wall jet that descends 
from the lid along the upper half of the downstream wall, with the periodic appearance 
and interaction of two small counter rotating tertiary vortices near the tip of this wall jet 
below where it errupts into the flow field as a whole. Note that the secondary vortex in the 
lower half of the cavity creates a relatively weak flow that rises along the lower half of the 
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downstream wall. The triggering instability that generates the asymptotic periodic flow 
is possibly the interaction of these two opposing flows. The relative strength of these two 
small tertiary vortices is related to the relative strength of the two opposing flows along 
the downstream wall. The second interesting feature of these dynamics is the evolution 
of the flow as the stronger upper small tertiary vortex is convected by the bounded shear 
layer across the middle of the cavity. The existence of local recirculation in the flow is 
possible because there are local streamfunction extremes in the center of the small vortices 
as they leave the wall, with saddle points in the streamfunction surface between the small 
vortices and the lower secondary vortex. The upper small tertiary vortex is convected to 
the streamfunction surface of the larger secondary vortex, the saddle point disappears, 
and the smaller vortex loses its independent identity to become a swelling on the side of 
the streamfunction surface of the larger vortex. This local swelling on the streamfunction 
surface is then convected across the bounded shear layer as a wavy disturbance in the 
streamfunction contour lines until the wall at the end of the shear layer is reached, where 
the local swelling is dissipated by the boundary layer. The bounded shear layer has an 
effect on the spectral signature of the local flow variables, since the nonlinear convection 
terms energize integer multiples of the basic frequency of the periodic shedding process at 
the end of the wall jet. In its final effects, the upper small vortex disturbs the flow all the 
way around the two large circulation patterns. Kinetic energy contours show this effect 
near the upper lid, and a close inspection of kinetic energy and vorticity surface plots shows 
an effect all along the upstream wall, and along both the lid and the lower wall. These 
effects away from the central band in the cavity are slight. The final interesting feature 
of these dyn am ics is the history of the lower small counterrotating tertiary vortex which 
originates a little later and slightly below the upper small vortex. This local flow pattern 
is convected away from the wall sooner than the upper small vortex, but it is weaker than 
the upper small vortex, and as the upper small vortex emerges into the flow the lower small 
vortex is pushed back toward the wall below its point of origin to be partially dissipated 
by the boundary layer. 

Figures lOa-i shows nine plots of data at the particular point with indexes (*,j) = 
(88,120), or (x, y) « (0.92,1.25), and for 4000 < t < 4100 . This point is near the end of 
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the wall jet that descends from the downstream end of the moving lid. The first three plots 
in Figures lOa-c are for the spectral densities of the point values of the streamfunction, 
the kinetic energy, and the vorticity. Figure lOd is a phase portrait of (^ n 31 ,V ; )> or 
the streamfunction point values plotted with a lag of 31 time steps. Figures lOe and 
lOf are similar time lagged phase portraits for the point values of the kinetic energy and 
the vorticity. Figures lOg-i are phase portraits of 2 {( tt ?,y) a + ( v *\y) 2 })’ 

where 10 is the vorticity, and (u”y,v"y), respectively. Figures lla-i shows nine similar plots 
at the particular point with indexes (t,j) = (32,106), or (x,y) « (0.33,1.10), and for 
4000 < t < 4100 . This point is approximately one third of the distance from the end of 
the bounded shear layer as it approaches the wall opposite the wall jet after crossing the 
center of the cavity. The data in Figures 10 and 11 is for every time step with a sample 
of 12800 time points for 4000 < t < 4100, with 295 iterations per cycle, and with slightly 
more than 43 cycles during this time interval. All of this data is represented in each of 
these plots. Notice the appearance of multiples of the fundamental harmonic in the data 
from the point near the end of the bounded shear layer. The point data in these figures 
shows the precisely periodic nature of the asymptotic flow with a fundamental frequency 
of approximately 0.433 « 

It was quite surprising to have to run the finer grid flow so long in order to obtain 
an asymptotic periodic state, since the coarse grid simulation converged at t ^ 1100. 
The recirculating boundary layer along the lower wall is present only in the refined grid 
asymptotic flow. The refined grid asymptotic solution has a larger total amount of kinetic 
energy on the interior grid than the coarse grid asymptotic flow . These differences between 
the two asymptotic flow solutions could be caused by a better resolution of gradients and 
momentum diffusion on the refined grid, particularly in the boundary layers. The precise 
long term repetition of the periodic asymptotic state, the general qualitative agreement 
between the asymptotic dynamics of the coarse and fine grid solutions, and the general 
qualitative agreement between the two coarse grid simulations at Re = 10000 with dif- 
ferent algorithms, all lead us to believe that the computed solution is a genuine periodic 
asymptotic state and not just a numerical artifact. 
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6: Remarks and Discussion 


We would like to include some general remarks about numerical results and method- 
ology as they pertain to this general new problem of computationally determining the 
qualitative nature of unsteady flows, especially for the long time behavior of unsteady 
flows. 

6.1: Methodology 

There are two views toward determining the qualitative behavior of a fluid flow. The 
first more predominant and easier view is to restrict attention to the branching diagram 
of the stationary steady flow problem. This is the bifurcation theory approach. The 
second view is to follow the unsteady flow to its final or asymptotic state. This is the 
dynamical systems approach. A main goal of both approaches is to determine all stable 
flow configurations. Both approaches may be said to be roughly equivalent if attention is 
restricted to stationary final states with u ( = 0. When progressing to the determination 
of qualitative flow behavior beyond the steady state the two approaches can yield different 
conclusions. 

In order to determine the transition to a periodic solution as the result of a Hopf 
bifurcation, the bifurcation theory approach typically uses an extended set of steady equa- 
tions by finding the eigenvalues of an associated Jacobian matrix, and then by finding 
the bifurcation parameter value at which a complex conjugate pair of such eigenvalues 
crosses the imaginary axis. Some arguments for the bifurcation theory steady equation 
analysis are: (l) computational cost is lower than an unsteady analysis; (2) the critical 
threshold parameter is predicted “exactly”, although after several mesh refinements; (3) 
the variation of the threshold parameter with other parameters can be studied in a similar 
way. 

Some arguments for using the dynamical systems approach of following an evolving 
flow to its asymptotic state are: (l) a time dependent computation “should” follow the 
“physics”, so stable branches can be explored with continuation procedures; (2) more 
complicated asymptotic states that occur further along a bifurcation diagram than the 
initial bifurcation points can be explored; (3) the dependence of the final state on initial 
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and intermediate conditions can be readily tested; (4) all of the intervening dynamics will 
be exhibited; (5) a complete picture of the final state is available within the limits of 
truncation error. We note that the unsteady dynamical systems approach as presently 
formulated does not readily pinpoint critical parameter threshold values for bifurcations. 

6.2: Solutions 

It is useful to draw some distinctions between various possible numerical flow solutions. 
A numerical algorithm for a fluid flow calculation is a system of algebraic equations, and 
this system of algebraic equations can have multiple solutions. Some of these solutions are 
valid and some are spurious. To quote Schreiber and Keller [18]: 

. . Furthermore, since the nonlinearity in the Navier-Stokes equations is quadratic, 
the approximating algebraic equations are also quadratic (in any reasonable scheme) . 
In the two-dimensional case with uniform mesh h in a domain of diameter 0(1) there 
are essentially N 2 = j-j unknowns and coupled quadratic equations. Now a basic 
result in algebraic geometry assures us that this algebraic system has 2^’ solutions, 
although some minor difficulties, i.e., “common intersection components,” must be 
eliminated or else there can be manifolds of solutions. If the flow problem of interest 
has a unique solution, we must hope that one of these 2 N solutions is a close approx- 
imation to it and that all of the others is spurious. This cursory account suggests that 
most of the numerical solutions are spurious! 

Fortunately most of the “numerical” solutions are also complex, so real computa- 
tions do not usually reveal them. Furthermore, solution procedures using continuation 
from known physical states may avoid them. But this is not always the case as we show 
in this note. Indeed even time marching schemes may lead to spurious steady states. 
Our results have revealed that this is particularly so when upstream differencing has 
been used in the driven cavity problem. 

Unfortunately there is at present no good theory to determine when a solution of the 
approximating problem is spurious and when it is “legitimate.” Indeed this imposes 
a severe burden on the computational fluid dynamicist to make additional tests on 
his results which will add weight to his assertion of their legitimacy. These tests 
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may affirm known physical or mathematical properties of the flow or else they may 
assure known approximation properties of the numerical method (i.e., h -truncation 
expansion, etc.).” 

We have addressed this issue by checking that our code and unsteady simulation approach 
accurately reproduce the known steady state solutions for the aspect ratio one driven cavity 
at Re = 5000 [6], and by checking that our simulations reproduce the same qualitative 
periodic flow dynamics at Re = 5000 on both coarse and refined grids. 

In a study of the aspect ratio one driven cavity, Glowinski, Keller, and Rinehart 

[5,p83l] have said: 

“A most interesting question is the possible occurence of multiple solutions as the 
Reynolds number increases beyond some critical value. . . . Actually and to our knowl- 
edge the computed solutions obtained in the range 0 < Re < 5000 by various authors 
using different methods agree quite well; this observation suggests that multiple solu- 
tions can only appear for greater values of Rt . 

The aspect ratio one driven cavity has corroborated solutions for 0 < Re < 10000 using 
both steady and unsteady algorithms. These various solutions differ from each other in 
terms of resolution and vortex details, but they are all steady solutions, and qualitatively 
similar in terms of the large scale vortex structures. 

An example of qualitative flow features that can distinguish solutions is given by 
Stokes flow with reasonable symmetry assumptions in lower cavity corners. A sequence 
of more than twenty corner vortices was resolved by Gustafson and Leben [10], where the 
smallest intensities were 0( lO' 80 ). This type of detail in a steady flow solution points 
out that numerical solutions can be distinguished by resolution which may not have any 
dynamical significance. There is a dynamically significant difference between steady and 
periodic asymptotic flow solutions. 

There are several known and mutually agreeing steady numerical solutions of the as- 
pect ratio two driven cavity for 0 < Re < 2000, plus the periodic solutions that we are 
reporting in this paper for Re = 5000. One of the features distinguishing the coarse grid 
and refined grid solutions that we have presented is the appearance in the refined mesh 
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solution of a recirculating boundary layer all across the bottom of the cavity. This sug- 
gests the possibility that a second grid dependent steady solution could appear before the 
critical point for the Hopf bifurcation. The grid dependent appearance of the recirculating 
boundary layer as a partially formed tertiary vortex in the bottom of the cavity does not 
seem to have a significant effect upon the essential periodic asymptotic dynamics in this 
Reynolds number range. 

6.3: Parameters 

Flows like the driven cavity can be thought of as three parameter bifurcation problems 
with a principal flow parameter such as Reynolds number, a principal geometry parameter 
such as aspect ratio, and a principal resolution or discretization parameter such as mesh 
size. Even though there are at least these three parameters that affect the cavity flow, we 
have concentrated here on the question of Hopf bifurcation with respect to the Reynolds 
number parameter with aspect ratio fixed at two. We conjecture that a Hopf bifurcation 
with respect to Reynolds number occurs in all lid driven cavities at all aspect ratios. If 
there is not a Hopf bifurcation for all aspect ratios, then it would be interesting to know 
the limits on aspect ratio for which there are Hopf bifurcations. If there are such limits, 
then we expect that they occur for small aspect ratios and not for large. A more general 
parametric investigation is underway. 
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7: Conclusions 


To summarize our numerical results for the aspect ratio two driven cavity. 

1] On a coarse 48 x 96 grid at Re = 10000 all measures indicate that the solution attains 
a periodic asymptotic state for t » 1500, with period 3.48 < T < 3.50. 

2] On a coarse 48 x 96 grid at Re = 5000 all measures indicate that the solution attains 
a periodic asymptotic state for t « 1100, with period 2.469 < T < 2.484. The 
asymptotic state for this flow is characterised by primary and secondary vortices in 
the upper and lower cavity, by periodic shedding of small tertiary counterrotating 
vortices off the downstream wall, by a wavy disturbance across the mid cavity shear 
layer, and by tertiary vortices in the bottom corners. 

3] On a refined 96 x 192 grid at Re = 5000 all measures indicate that the solution 
attains a periodic asymptotic state for t » 3700, with period 2.305 < T < 2.313. 
The asymptotic state for this flow is characterised by primary and secondary vortices 
in the upper and lower cavity, by periodic shedding of small tertiary counterrotating 
vortices off the downstream wall, by a wavy disturbance across the mid cavity shear 
layer, by a partially formed third principle vortex at the bottom of the cavity, and by 
quartiary vortices in the bottom corners. 

Based upon these computations we conjecture that the Navier-Stokes equations for the 
aspect ratio two driven cavity possess a Hopf bifurcation in the interval 2000 < Re < 5000, 
since the transition has been computationally demonstrated. 

Many interesting flow dynamics are shown by simulating the full time dependent flow 
history from the initial no flow state to the final periodic state shows. Particular examples 
are the development from early time of the periodic shedding of counterrotating vortex 
couples from the downstream wall, and the two unexpected transient oscillatory regimes 
before converging to a permanent periodic solution in the Re = 5000 calculation with a 
96 x 192 grid. Following the dynamical history of the flow helps in the interpretation of 
the final asymptotic state. 

The investigation of the qualitative properties of unsteady flows treated as infinite 
dimensional dissipative dynamical systems is ushering in a new chapter of numerical anal- 
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ysis. We have presented a number of qualitative measures and indicators that we have 
found to be useful in this type of study. In particular, we stress the dangers of relying 
totally upon any one of these measures. 
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